## Copyright (C) 2012 Estela Moura
## 
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## candidato_ACO-TSP
##Retorna: BestTourLength
##Recebe: distances_matrix, number_of_ants, MaxIterations, par_alpha, par_beta, rho, target_length

## Author: Estela Moura <estela@kubuntu-Inspiron>
## Created: 2012-06-13

function [BestTourLength, BestTour] = ACO_TCP_from_Internet_adaptado(cities_matrix, number_of_ants, MaxIterations, par_alpha, par_beta, rho)

% Ant Colony System for the TSP
%
% ACSTSP(distances_matrix, number_of_ants, MaxIterations, par_alpha,par_beta, rho, target_length)
%
% distances_matrix: symmeric real (NrOfNodes x NrOfNodes)-matrix containing distances
% between all nodes: d[i,i] = 0, d[i,j]=d[j,i]
% MaxIterations: maximum number of iterations to run
%
% parameter standard values:
% number_of_ants = NrOfNodes
% par_alpha = 1
% par_beta = 9
% rho = 0.9
%
% returns [BestTourLength]

for(i=1:size(cities_matrix, 1))
    for(j=1:size(cities_matrix, 1))
      if(i==j)
	distances_matrix(i,j)=0;
      else
	distances_matrix(i,j)=sqrt((cities_matrix(j,2)-cities_matrix(i,2))^2 + cities_matrix(j,3)-cities_matrix(i,3)^2);
      endif
    endfor
  endfor

d = distances_matrix; %define d
n = max(size(d));

m = number_of_ants;
t_max = MaxIterations;

[L_nn, P_nn] = NearestNeighborTSP(d);

L_best = inf;
T_best = 0;

% INITIALIZATION ===========================================================

% pheromone trails
c = 1 / (n * L_nn);
tau = ones(n,n) * c;

% place m ants in n nodes
ant_tours = zeros(m, n+1);
# # # ant_tours(:,1) = randint(m,1,[1,25]); -> esta função "randint" não existe no Octave
# # # 	aqui, cria-se uma matriz mx1 com valores entre 1 e 25, uma solução equivalente em Octave está abaixo:
ant_tours(:,1)=randi(25, m, 1);

t = 1;
while (t <= t_max)

% CREATE TOURS =============================================================

	for s = 2 : n
		for k = 1 : m
            current_node = ant_tours(k,s-1);
            visited = ant_tours(k,:);
            to_visit = setdiff([1:n],visited);
            c_tv = length(to_visit);
            p = zeros(1,c_tv);
            for i = 1 : c_tv
                p(i) = (tau(current_node,to_visit(i)))^par_alpha * (1/d(current_node,to_visit(i)))^par_beta;
            end
            sum_p = sum(p);
            p = p / sum_p;
            for i = 2 : c_tv
                p(i) = p(i) + p(i-1);
            end
            r = rand;
            select = to_visit(c_tv);
			for i = 1 : c_tv
                if (r <= p(i))
                    select = to_visit(i);
                    break;
                end
			end
            city_to_visit = select;
            ant_tours(k,s) = city_to_visit;
            tau(current_node,city_to_visit) = (1 - rho) * tau(current_node,city_to_visit) + c;
		end
	end

% UPDATE ===================================================================

	ant_tours(:,n+1) = ant_tours(:,1);
	L_T = zeros(1,m);
	best_ant = 1;
	for k = 1 : m
        P = ant_tours(k,:);
        L = 0;
        for i = 1 : n
            L = L + d(P(i),P(i+1));
        end
        L_T(k) = L;
        if (L_T(k) < L_T(best_ant))
            best_ant = k;
        end
	end
	L_min = min(L_T);
	T_min = ant_tours(best_ant,:);

% update pheromone trails;
	for i = 1 : n
        tau(T_min(i),T_min(i+1)) = (1 - rho) * tau(T_min(i),T_min(i+1)) + rho / L_min;
	end

% COMPLETE =================================================================
	clc;
	t = t + 1;
	current_cities = ant_tours(:,n);
	ant_tours = zeros(m, n+1);
	ant_tours(:,1) = current_cities;
	if (L_min < L_best)
        L_best = L_min;
        T_best = T_min;
	end

end % ends while
BestTourLength=L_best;
BestTour = T_best;

endfunction
